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Abstract 

This document describes a grid generation program for use with the CAP-TSD 
transonic small disturbance code* The program runs interactively in FORTRAN on the 
Sim Workstation. A fifth-degree polynomial is used to map the grid index onto the 
computational coordinate. The grid is plotted to aid in the assessment of its quality and 
may be saved on a file in NAMELIST format. 

Approach 

In the transonic small disturbance code CAP-TSD 1 the computational grid consists 
of independent meshes in each of the computational coordinates XI, ETA, and ZETA with 
indices I, J, and K. That is, XI(I) is a one- dimensional array which maps the integers 
I = 1 to NXT onto the computational coordinate XI, with similar expressions for the 
other two arrays. Typically, the streamwise XI and vertical ZETA meshes extend 10 to 
40 chord lengths from the configuration and the spanwise ETA mesh extends 2 to 4 span 
lengths. 

Each of the meshes may be built up from a sequence of independent meshes for each 
region of the configuration. For a simple wing, the XI mesh might consist of three regions: 
(1) a mesh extending from the upstream boundary of the computational domain to the 
wing leading edge; (2) one from the leading edge to the trailing edge; (3) one from the 
trailing edge to the downstream boundary. It is desirable that configuration changes such 
as wing, tail, and body edges lie either on or half-way between the mesh points. For the 
simple wing, the up and downstream far-field boundaries would lie on mesh points with 
indices 1 and NXT, the wing leading edge would lie half-way between two mesh points, 
and the trailing edge would coincide with a mesli point. 
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In general, the mesh spacing should be finer at. region boundaries and in other regions 
of rapid changes in the flow. Typically, mesh spacing on a lifting surface might range 
from 0.2 to 2.0 percent chord in XI and from 0.5 to 5.0 percent chord in ETA. The 
ZETA mesh spacing at the lifting surface should be comparable to that used cliordwise 
with the addition of two points lying very close above and below the surface at 0.001 
percent chord. The non-reflecting boundary conditions will be better approximated by 
the difference equations if the mesh spacing is kept below one chord length at all far-field 
boundaries. 

Method 

In this grid generation program, the mesh for each independent region is constructed 
from a fifth-degree polynomial which furnishes the computational coordinate (called X in 
the program) as a function of the index of the mesh point (called I). This quintic polyno- 
mial has six independent constants which are chosen from the following three conditions 
at each endpoint of the region: (1) the mesh index at the computational coordinate of 
the region boundary; (2) the mesh spacing at the boundary; (3) zero second derivative 
of the polynomial at the boundary. This last condition insures mesh smoothness at the 
boundaries. The index at the region boundary would be an integer for the case in which 
the mesh point lies on the boundary (e.g., at the wing trailing edge), but would be a 
half-integer for the case in which the boundary falls between mesh points (e.g., at the 
leading edge). 

Although the meshes for each region may be computed independently, one should 
maintain continuity and smoothness at the region boundaries by using the same mesh 
spacing at a common boundary. The program facilitates this condition by allowing for 
the computation of the mesh for more than one region at a time. 

Computer Program 

The computer program was originally written in BASIC for the IBM-PC with EGA 
graphics card. The program described herein is written in FORTRAN for the Sun Work- 
station using Sun CGI graphics. A fisting of the program is given in the Appendix. 

The program begins by displaying a screen of general information. You are asked to 
enter (Return) to continue. The program then asks for the number of regions in the mesh 
you wish to design. As with all the numeric input, a default is provided. You are next 
asked for three quantities at each region boundary: (1) the computational coordinate of 
the boundary; (2) the mesh spacing at. the boundary; (3) the mesh index of the boundary 
point (integer or half- integer, greater than zero). In each case, the current value of the 
quantity is displayed and may be left unchanged by entering (Return). 

When the entry of the data for each boundary is complete, the program writes 
these defining values on the file “param.teni”. The program then fists the mesh index, 
coordinate, first, and second differences on the screen. Negative values for the first 
difference indicate that the mesh has folded on itself — fewer points in the region or 
smaller spacings at the boundaries will fix this. 


2 


Response to the next prompt will clear the screen mid plot the mesh as horizontal 
lines against X. The mesh index I is shown on the horizontal axis and the mesh func- 
tion, first difference, and second difference are plotted against the index. The defining 
parameters for the mesh and its maximum and minimum spacings are printed at the top 
of the screen. The region boundaries are shown as tic marks with the boundary numbers 
along the right edge of the frame. The plot is in color; a black-on-wliite version may 
be obtained by changing “call colortable8” to “call colortable2” in the source code. To 
clear the plot screen you must enter “go” and (Return). Any other response will exit the 
program. 

After exiting from the plot, you are asked whether you wish to save this grid on a 
file. This option will write the mesh coordinates in a five decimal place format, five to 
a line separated by commas, on your named file. Later, this file may be edited into the 
CAP-TSD input namelist. You are asked next if you wish to exit the program. If your 
mesh was not satisfactory, you may respond to the prompt simply with (Return) and the 
program will re-enter its input phase. Your current values will be displayed and you may 
make any desired changes. 

Output of Computer Program 

As examples of the mesh generation program output, results for the five default grids 
are shown. The meshes are described in terms of applications to particular configurations. 
However, these meshes are too coarse and the far-field boundaries are too close to provide 
accurate results. 

For each example, a screen dump of the black-on-white plot of the mesh is shown. 
The command “screendump | rasfilter8tol | lpr -s -v typed in a window other than 
the one containing the plot, will dump the entire screen to the printer. The screen used 
in the examples is one with the largest square area for plotting available on the Sun 
Workstation. It is 900 pixels tall and 880 pixels wide. Use of too small a window will 
lead to distortions in the lettering and other problems. 
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The default mesh for one region is shown in figure 1. This mesh is intended to 
illustrate vertical ZETA mesh properties for an isolated wing. The spacing is fine near 
the wing surface (X = 0) and stretches smoothly to a large spacing at the boundary 
(X = 8). The grid function X is plotted with long dashes, the spacing dX with short 
dashes, and the second difference ddX with dot-daslies. The origin for the latter curve is 
arbitrary, but this difference does approach zero at the region boundaries. For accurate 
calculations with CAP-TSD, the mesh would need two to three times the number of 
points and mesh extent as illustrated here. 



Figure 1.- Plot of default mesh with 1 region. 
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The default mesli for two regions is shown in figure 2. This mesh is typical of a very 
coarse spanwise ETA mesh for a simple wing. The spacing decreases smoothly from the 
root (X = 0) to the tip (X = 1) and then increases to the far-field boundary (X = 3). 
The wing tip lies between two mesh points (X = 10.5). The number of points on the 
wing was chosen to provide nearly uniform spacing near the root. The number of points 
off the wing was chosen to provide similar stretchings about, the tip (the curve for dX is 
symmetrical there) and to achieve the maximum spacing at the far-field boundary. All 
of the curves possess the nice property of smoothness. 



Figure 2- Plot of default mesh with 2 regions. 
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The default mesh for three regions is shown in figure 3. Tliis mesh illustrates the 
features of a streamwise XI mesh for an isolated wing, but it has too few points, too close 
far-field boundaries, and too large a spacing on the wing to be practical. The wing lies 
in the middle region with leading edge at X = 0 and trailing edge at X = 1. The leading 
edge lies between mesh points 13 and 14 and the trailing edge at point 28. The mesh is 
stretched somewhat over the middle of the wing chord and has its finest spacing at the 
wing edges. It achieves maximum spacing at the up and downstream boundaries. 



CAPgrid for 3 regions Itot = 41 dXmin = 0.05000 dXmax = 0.40000 

Xbd = -3.00000 0.00000 1.00000 4.00000 X 

dXbd = 0.40000 0.05000 0.05000 0.40000 dX 

Ibd = 1.0 13.5 28.0 41.0 ddX 
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Figure 3.- Plot of default mesh with 3 regions. 
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The default mesli for four regions is shown in figure 4. This example illustrates a 
coarse spanwise ETA mesh for a fuselage/wing/tail configuration. The first mesh point 
lies on the centerline at X = 0. The side of the fuselage is at X = 0.1 and lies between 
mesh points. The tail tip is at X = 0.5 and the wing tip at X = 1.0, with both lying 
between mesh points. The far-field boundary is at X = 2.0. The mesh is clustered 
somewhat near the lifting surface tips with spacings of three percent of the configuration 
semispan, and stretches symmetrically away from them. There are three mesh points on 
the (half) fuselage, lxine on the tail, and 21 on the wing. 



CAPgrld for 4 regions Itot = 34 dXnnn = 0.03000 dXmax = 0.20000 

Xbd = 0.00000 0.10000 0.50000 1.00000 2.00000 X 

dXbd = 0.04000 0.04000 0.03000 0.03000 0.20000 dX 

Ibd - 1.0 3.5 12.5 24.5 34.0 ddX 
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Figure 4.- Plot of default mesh with 4 regions. 
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Conclusion 


An interactive program to aid in designing grids for use with the CAP-TSD com- 
puter code has been described. The program uses a fifth-degree polynomial to define the 
computational coordinate as a function of the grid point index. The program default 
grids, which illustrate some of the considerations in producing suitable grids for use with 
CAP-TSD, are shown. 
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APPENDIX - Program Listing 


program capgrid 
c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c C 
c Grid generation for multiple regions with restart c 
c c 
c SUN FORTRAN version - 90/08/09 c 
c c 
c The compile statement is: c 
c c 
c f77 -o capgrid capgrid. f -lcgi77 -lcgi -lsunwindow -lpixrect -1m c 
c c 
c Samuel R. Bland UAB/SDyD (804) 864-2272 c 
c NASA Langley Research Center c 
c c 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

common / bound / nreg, xb(0:5), rib(0:5), db(0:5), kb(0:5), 

$ rir(0:5) / irm(0:5), ril(0:5), ilp(0:5), 

$ xr (0:5) , xrm(0:5), xl(0:5), xlp(0:5) 

common / mesh / x(999), dx(999), ddx(999), 

$ imin, imax, itot, xmin, xmax, dxmin, dxmax 


c 


c 

c 


c 


c 


c 

c 

c 

c 

c 

c 

100 

c 

c 

c 


character zzz$*8 
character yn$*8 
character gridf ile$ *40 

integer kb(0:5), ilp(0:5) r irm(0:5) 

real ril(0:5) / rir(0:5) 

real xb(0:5), rib (0:5), db(0:5) 

real xl(0:5), xlp(0:5), xrm(0:5), xr(0:5) 

real*8 x(999), dx(999), ddx(999) 
real*8 c ( 6 ) 

jmax - 6 
nreg — 2 
iter = 0 

call descrip 

call colortable8 
call colortable2 

Begin grid iteration 

continue 

call enterinput 

call mapping (jmax,c) 

call diffs 
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c 

c 

c 


call printscrn 

call plotgrid 

print ' ( // , a , $ ) ' , 
read ' (a8 ) ' , yn$ 


' Save the grid on a fide (y/n)? 


if (yn$ .eq. ' y' ) then 
print ' (a) ' , 


Input the file name' 


c 

c 


read '(a40)', gridfile$ 
call writegrid (gridfile$) 


print ' (/ ) ' 

print ' (a,$)', ' Enter 
read ' (a8 ) ' , zzz$ 
endif 


<Return> to continue 


print ' (//,a,$)', ' Do you want to EXIT the program (y/n)? ' 

read ' (a8) ' , yn$ 

if (yn$ .ne. 'y') then 
goto 100 
endif 


END OF PROGRAM 




print ' (/) ' 

P ^ ^ ^ **★★★★****★*★* 

print ' (/ ) ' 

c 

end 

subroutine enter input 
c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 


Enter the region boundary coordinate, grid spacing, and 

boundary index. Compute the grid coordinate and index for 
each region . 

The definitions are: 

Region boundary index (0-nreg) 

Coordinate of the region boundary 
Grid spacing at the region boundary 
Real grid index at region boundary (integer or 
half integer 

Flag to indicate whether boundary point is a 
grid point 

Real index of left boundary point 
Real index of right boundary point 
Index of leftmost point wholly within region 
Index of rightmost point wholly within region 

Coordinate of left boundary point 
Coordinate of right boundary point 

Coordinate of leftmost point wholly within region 
Coordinate of rightmost point wholly within regi 


c 

n 

c 

xb (n) 

c 

db (n) 

c 

rib (n) 

c 


c 

kb (n) 

c 


c 


c 

ril(n) 

c 

rir (n) 

c 

ilp(n) 

c 

irm(n) 

c 


c 

xl(n) 

c 

xr (n) 

c 

xlp(n) 

c 

xrm(n) 


c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

on c 
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C c 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 


c 


c 


$ 

$ 

$ 


common / bound / 
common /mesh / 


nreg, xb(0:5), rib(0:5), db(0:5), kb(0:5), 
rir(0:5), irm(0:5), ril(0:5), ilp(0:5), 
xr (0 : 5) , xrm(0:5), xl(0:5), xlp(0:5) 
x ( 999) , dx (999) , ddx(999), 

imin, imax, itot, xmin, xmax, dxmin, dxmax 


character*8 inp$ 
character*l n$ 

real*8 x(999), dx(999), ddx(999) 


iter = iter + 1 


c 


c 
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c 


print ' (/) ' 
print ' (a) ' , 
print ' (/) ' 
print ' (a) ' , 
print ' (a) ' , 
print ' (/) ' 


************ CAP-TSD GRID DESIGN 


★ ■*■■*■★**•*•*•★* ★ k t 


You may leave 
responding 


any parameter unchanged by' 

to the input request with <Return>' 


print ' (a,i3,a,$)', ' Number of regions =' , nreg, ' ? ' 
read ' (a8 ) ' , inp$ 
iflag = 0 
do 10 k = 0, 9 

if (iflag .eq. 0) then 
write (n$, ' (il ) ' ) k 
if (index (inp$, n$) .ne. 0) then 
iflag = 1 
iter = 1 

read (inp$, '(i8)') nreg 
endif 
endif 
continue 


if (iter .eq. 1) then 
call iniparam 
endif 


c 

do 400 n = 0, nreg 

print '(/,a,i2)', ' Data at boundary n 

c 

print ' (a,fl0.5,a,$) ' Coordinate =' , xb(n), ' ? ' 

read ' (a8 ) ' , inp$ 
if (index (inp$ ) .ne. 0) then 
read (inp$, '(fl0.5)') xb(n) 
else 

iflag = 0 
do 100 k = 0, 9 

if (iflag .eq. 0) then 
write (n$, ' (il) ' ) k 
if (index (inp$, n$ ) .ne. 0) then 
iflag = 1 

read (inp$, '(i8)') iin 
xb(n) = iin 
endif 
endif 

100 continue 

endif 
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c 

print ' (a,fl0.5,a,$) ' , ' Spacing =' , db(n), ' ? ' 

read ' (a8 ) ' , inp$ 

if (index (inp$, '.' ) .ne. 0) then 
read (inp$, '(fl0.5)') db(n) 
else 

iflag = 0 
do 200 k = 0, 9 

if (iflag .eq. 0) then 
write (n$, ' (il) ' ) k 
if (index (inp$, n$) .ne. 0) then 
iflag = 1 

read (inp$, '(i8)') iin 
db(n) = iin 
endif 
endif 

200 continue 

endif 
c 

print ' (a,fl0.1,a,$)', ' Index = ' , rib(n), ' ? ' 

read ' (a8) ' , inp$ 
if (index (inp$, '.' ) .ne. 0) then 
read (inp$, '(fl 0 . 5 )') rib(n) 
else 

iflag = 0 
do 300 k = 0, 9 

if (iflag .eq. 0) then 
write (n$, ' (il) ' ) k 
if (index <inp$, n$) .ne. 0) then 
iflag = 1 

read (inp$, '(i8)') ii n 
rib(n) = iin 
endif 
endif 

300 continue 

endif 

c 

if (rib (n) .eq. int(rib(n))) then 
kb(n) = 1 
else 

kb (n) = 2 

rib(n) = int(rib(n)) + .5 
endif 

c 

400 continue 
c 

if (kb(0) .eq. 1) then 
imin = rib(0) 
else 

imin = rib(0) + .5 
endif 
c 

if (kb(nreg) .eq. 1) then 
imax = rib(nreg) 
else 

imax = rib(nreg) - .5 
endif 
c 
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c 

c 

c 

c 

c 

c 

c 


500 

c 


itot * imax - imin + 1 
xmin = xb(0) 
xmax = xb(nreg) 

Save the parameters that define the grid on file "par am . t < •in" 


call outparam 

print '(//,a)' f ' Grid definition is saved on file "param. tern" ' 

Compute region boundary parameters 


do 500 n - 0, nreg 
rir(n) = rib(n) 
ril (n + 1) = rib (n) 
xr(n) - xb(n) 
xl (n + 1) =* xb (n) 
if (kb (n) .eq. D then 
irm(n) * rib(n) - 1. 


ilp (n+1) 
xrm(n) = 
xlp (n+1 ) 
else 

irm(n) = 
ilp (n+1) 
xrm(n) = 
xlp (n+1 ) 
endif 
continue 


rib(n) + 1. 
xb ( n ) - db ( n ) 

- xb(n) + db(n) 

rib(n) - .5 
= rib (n) + .5 
xb(n) - .5 * db(n) 

= xb ( n ) + .5 * db ( n ) 


return 

end 


subroutine diffs 
c 

ccccccccccccccccccccccccccccccccccccccccccccccc 

c c 

c Compute the 1st and 2nd differences c • 

c c 

ccccccccccccccccccccccccccccccccccccccccccccccc 

common / mesh / x(999), dx(999), ddx(999), 

$ imin, imax, itot, xmin, xmax, dxmin, dxmax 

c 

real*8 x(999), dx(999), ddx(999) 

c 

do 100 i * imin, imax - 1 
dx (i) = x ( i+1 ) - x (i ) 

100 continue 
c 

do 110 i = imin + 1, imax - 1 

ddx(i) = x (i+1 ) - 2. * x(i) + x(i-l) 

110 continue 
c 


14 



120 


dxmin = 1000. 
dxmax = 0 . 

do 120 i = imin, imax - 1 

if (dx(i) .It. dxmin) then 
dxmin - dx(i) 
endif 

if (dx(i) .gt. dxmax) then 
dxmax = dx(i) 
endif 
continue 


return 

end 

subroutine printscrn 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccc 


c 

c 

c 


Print the mesh values and differences 


c 

c 

c 


ccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

common / mesh / x(999), dx(999), ddx(999), 

$ imin, imax, itot, xmin, xmax, 


c 

c 

c 


real*8 x(999), dx(999), ddx<999) 
print ' (/, a) ' , ' I X (I) 


dX ( I ) 


100 


print 9 (i3, 2f 12 . 5 ) 9 , imin, x(imin), dx(imin) 
do 100 i = imin + 1, imax - 1 

print 9 ( i3, 3f 12 . 5 ) 9 , i, x(i), dx(i), ddx(i) 
continue 

print r (i3, f 12 . 5, /) 9 , imax, x(imax) 

return 

end 


subroutine colortable2 


ccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 
c Define a 2 color table with the values: c 
c c 
c 0 white 1-7 black c 
c c 
ccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

common / color / lred, lgrn, lblu, ncolors 


c 

c 


integer lred(8), lgrn(8), lblu(8) 
ncolors * 8 


lred ( 1) = 255 
lgrn ( 1) = 255 
lblu (1) = 255 
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dxmin, dxmax 
ddX ( I ) ' 



100 

c 


do 100 i = 2, 8 
lred(i) = 0 
lgrn(i) = 0 
lblu(i) = 0 
continue 

return 

end 


subroutine colortable8 
c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c 





c 

c 

Define . 

an 8 color 

table 

with the values: 

c 

c 





c 

c 

0 

dark blue 

4 

red 

c 

c 

1 

blue 

5 

magenta 

c 

c 

2 

green 

6 

yellow 

c 

c 

3 

cyan 

7 

white 

c 

c 





c 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

common / color / Ired, lgrn, Iblu, ncolors 
c 

integer lred(8), lgrn(8), lblu(8) 
integer lc(2) 
c 

ncolors - 8 
lc(l) = 0 

lc (2) = 255 
c 

i = 0 

do 300 ir ■ 1 # 2 
do 200 ig * 1, 2 
do 100 ib = 1, 2 
i = i + 1 
lred(i) - lc(ir) 
lgrn(i) =* lc(ig) 
lblu(i) = lc (ib) 

100 continue 

200 continue 

300 continue 
c 

lblu(l) = 160 

c 

return 

end 


subroutine plot grid 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c C 
c Open the CGI graphics and plot the grid, the grid function, c 
c and the 1st and 2nd differences. c 
c c 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 
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c 


c 


c 

c 


c 

c 

c 

c 

c 

c 


100 

c 

c 

c 


110 


common / bound / nreg, xb(0:5), rib (0:5), db(0:5), kb (0:5), 

rir (0:5) , irm(0:5), ril(0:5), ilp(0:5), 
xr (0 : 5) , xrm(0:5), xl(0:5), xlp(0:5) 
common / mesh / x(999), dx(999), ddx(999), 

imin, imax, itot, xmin, xmax, dxmin, dxmax 
common / color / lred, lgrn, lblu, ncolors 
common / plot / lxof, lyof, rxs, ris, rxds, rxdds 

character*80 ntex$ 
character*8 zzz$ 

integer lred(8), lgrn (8), lblu (8) 
integer lxp(999), lyp(999) 
real*8 x(999), dx(999), ddx(999) 


print ' (a)', ' The program will -plot the grid next. To end the' 

print ' (a)', ' plot without exiting the program type "go"' 

print ' (a) ' , ' then <Return> after the plot is completed.' 

' (/ r $) ' r ' Enter <Return> to continue ' 
read ' (a8 ) ' , zzz$ 
print ' (/) ' 


call 

call 

call 

call 

call 

Draw 


cfopencgi () 

cf openvws (name, 0,0,0, 1 , CGPIXWINDD, ncolors, 
'Color', 0,0, 0,0, 5,0) 

cf cotable (0, lred, lgrn, lblu, ncolors) 

cf lnspecmode (0) 

setupplot 

the grid 


call cflnwidth (1.) 
call cflntype (0) 
call cflncolor (2) 
do 100 i ~ imin, imax 
lxp(l) = lxof 

lxp ( 2 ) - lxp(l) + ris * (itot - 1) 
lyp(l) = lyof + rxs * (x(i) - xmin) 
lyp (2) = lyp(l) 
call cfpolyline (lxp, lyp, 2) 
continue 


Plot X 


call cflnwidth (2.) 
call cflntype (5) 
call cflncolor (7) 
ip = 0 

do 110 i ~ imin, imax 
ip = ip + 1 

lxp ( ip) = lxof + ris * (i - imin) 
lyp(ip) = lyof + rxs * (x(i) - xmin) 
continue 

call cfpolyline (lxp, lyp, ip) 
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Plot dx 


c 

c 

c 


120 


c 

c 

c 


130 


c 

c 

c 


200 


210 

c 

c 

c 

c 


call cf Inwidth (2.) 
call cflntype (2) 
call cflncolor (3) 
ip = 0 

do 120 i = imin, imax - 1 
ip = ip + 1 

lxp(ip) = lxof + ris * (i — imin + .5) 
lyp(ip) = lyof + rxds * dx(i) 
continue 

call cfpolyline (lxp,lyp,ip) 

Plot ddX 

call cflnwidth (2.) 
call cflntype (3) 
call cflncolor (6) 
ip = 0 

do 130 i = imin + 1, imax - 1 
ip = ip + 1 

lxp(ip) = lxof + ris * (i - imin) 

lyp(ip) = lyof + rxs * ( (xmax - xmin) * .75) + rxdds * ddx(i) 
continue 

call cfpolyline (lxp,lyp,ip) 

Draw and label region boundaries 

call cflnwidth (1.) 
call cflntype (0) 
call cflncolor (5) 
call cftextcolor (7) 
call cftextfontix (4) 
lxt = 32000 
do 210 n = 0, nreg 

lyp(l) = lyof + rxs * (xb(n) - xmin) 

lyp (2) = lypd) 

do 200 nl = 0, nreg 

lxp(l) = lxof + ris * ( rib (nl ) - imin) - 400 
lxp(2) = lxp(l) + 800 
call cfpolyline (lxp,lyp,2) 
continue 

lyt = lyp(l) - 180 

write (ntex$. Mil) ' ) n 

call cftext (lxt, lyt, ntex$) 

lxp(l) = lxof + ris * (rib (n) - imin) 

lxp(2) = lxp(l) 

lyp ( 1 ) = lyof + rxs * (xb(0) - xmin) 
lyp (2) = lyof + rxs * (xb(nreg) - xmin) 
call cfpolyline (lxp,lyp,2) 
continue 

pause 

call cfclosevws (name) 
call cfclosecgi () 


return 

end 
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subroutine setupplot 


c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c c 

c Define the plot origin and scale factors and draw and label c 

c the axes . c 

c c 

c The plot units are in the range 0-32767 for x and y c 

c and use the largest square available in the current window. c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c 


c 


c 


c 

c 

c 


c 

c 

c 


c 


c 


common / bound / 

$ 

$ 

common / mesh / 

$ 

common / color / 
common / plot / 


nr eg, xb(0:5), rib (0:5), db(0:5), kb (0:5), 
rir(0:5), irm(0:5), ril(0:5), ilp(0:5), 
xr (0 : 5) , xrm(0:5), xl(0:5), xlp(0:5) 
x (999) , dx (999) , ddx(999), 

imin, imax, itot, xmin, xmax, dxmin, dxmax 

lred, lgrn, lblu, ncolors 

lxof, lyof, rxs, ris, rxds, rxdds 


character*80 ntex$ 

integer lred(8), lgrn(8), lblu(8) 
integer lxp(999), lyp(999) 
real*8 x(999), dx(999), ddx(999) 


lxof = 5200 
lyof = 2300 

rxs = 27000. / (xmax - xmin) 
ris = 25900. / (imax - imin) 
rxds = 4 . * rxs 
rxdds = 8 . * rxs 

Draw rectangle around plotting area 

call cfperimspecmode (0) 
call cfperimwidth (0.) 
call cfperimtype (0) 
call cfperimcolor (4) 

call cf rectangle (32717,50,50,32717) 

Write legend on plot 

call cftextcolor (7) 
call cftextfontix (4) 


lxt = 800 
lyt = 31900 

write (ntex$, ' (a, il, a, i4 , 3x, a, f 7 . 5, 3x, a, f 8 . 5 ) ' ) 

$ ' CAPgrid for ', nreg, ' regions Itot =' , 

$ itot, ' dXmin « ', dxmin, 'dXmax = ' , dxmax 

call cf text (lxt, lyt, ntex$) 


lxt = 1000 
lyt = lyt - 800 

call cftext (lxt, lyt,' Xbd =' ) 


c 


lyt = lyt - 600 

call cftext (lxt, lyt, ' dXbd =' ) 



c 

lyt = lyt - 600 
call cftext (lxt, lyt,' Ibd =' ) 
c 

lxt = lxt + 2000 
lyt = lyt + 2 * 600 

write (ntex$, ' (6 (f 11 . 5) ) ' ) (xb(n), n = 0, nreg) 

call cftext ( lxt, lyt , ntex$ ) 
c 

lyt = lyt - 600 

write (ntex$, ' (6 (f 11 . 5) ) ' ) (db(n), n = 0, nreg) 

call cftext (lxt, lyt, ntex$) 
c 

lyt = lyt - 600 

write (ntex$, ' (6 (f 7 . 1, 4x) ) ' ) (rib(n), n = 0, nreg) 
call cftext (lxt, lyt, ntex$) 
c 

c Draw key on plot 

c 

lxt = 30400 
call cflnwidth (2.) 
lxp(l) = 27500 
lxp (2 ) = 30000 

c 

lyt = lyt + 1200 

call cftext (lxt, lyt,' X') 

lyp(l) = lyt + 180 

lyp (2) - lyp(l) 
call cflntype (5) 
call cflncolor (7) 
call cfpolyline (lxp, lyp, 2) 
c 

lyt = lyt - 600 
call cftext (lxt, lyt,' dX' ) 
lyp (1) = lyt + 180 
lyp (2) = lyp(l) 
call cflntype (2) 
call cflncolor (3) 
call cfpolyline (lxp, lyp, 2) 
c 

lyt = lyt - 600 
call cftext (lxt, lyt, ' ddX' ) 
lyp(l) = lyt + 180 
1 yp ( 2 ) = lyp (1) 
call cflntype (3) 
call cflncolor (6) 
call cfpolyline (lxp, lyp, 2) 
c 

c Draw the I-axis 

c 

call cflnwidth (1.) 
call cflntype (0) 
call cflncolor (7) 
lxp(l) = lxof 

lxp (2) = lxof + ris * (imax - imin) 

lyp(l) = lyof - 1000 

lyp (2 ) = lyp(l) 

call cfpolyline (lxp, lyp, 2) 
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c 


do 200 i = imin, imax 

lxp(l) = lxof + ris * (i - imin) 

lxp { 2 ) = lxp(l) 

lyp(l) = lyof - 1000 

lyp (2) = lyp(l) + 500 

call cfpolyline (lxp, lyp, 2) 

200 continue 
c 

c Label the I-axis 

c 

call cftextcolor (7) 
call cftextfontix (4) 
lyt = 400 

lxt = lxof + ris * .5 * (imax - imin) - 100 
call cftext (lxt, lyt, ' I' ) 
c 

lxt = lxof - 700 
lyt = 600 

write (ntex$, ' (i3) ' ) imin 
call cftext (lxt, lyt , ntex$ ) 
c 

lxt = lxof + ris * (imax - imin) - 700 
write (ntex$, ' (i3) ' ) imax 
call cftext (lxt, lyt, ntex$) 

c 

c Draw the X-axis 

c 

call cflnwidth (1.) 
call cflntype (0) 
call cflncolor (7) 
lxp(l) = lxof - 1500 
lxp (2) = lxp ( 1 ) 
lyp(l) = lyof 

lyp (2) = lyp(l) + rxs * (xmax - xmin) 
call cfpolyline (lxp, lyp, 2) 

c 

do 300 j = 0, 1 

lxp { 1 ) = lxof - 1500 
lxp (2) = lxp ( 1 ) + 500 

lyp (1) — lyof + j * rxs * (xmax - xmin) 

lyp (2) = lyp ( 1 ) 

call cfpolyline (lxp, lyp, 2) 

300 continue 
c 

do 310 xt = - 100., 100. 

if ( (xt .It. xmax) .and. (xt . gt . xmin)) then 
lxp ( 1 ) = lxof - 1500 
lxp (2 ) = lxp ( 1 ) + 500 
lyp(l) = lyof + rxs * (xt - xmin) 
lyp (2) = lyp (1) 
call cfpolyline (lxp, lyp, 2) 
endif 
310 continue 
c 
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c Label the X-axis 

c 

call cftextcolor (7) 
call cftextfontix (4) 
lxt — 400 
c 

lyt “ lyof + rxs * .5 * (xmax - xmin) - 200 
call cftext (lxt, lyt,' X') 

c 

lyt = lyof - 200 
write (ntex$, ' (f 9 . 4 ) ' ) xmin 
call cftext (lxt, lyt, ntex$) 
c 

lyt * lyof + rxs * (xmax - xmin) - 200 
write (ntex$, ' (f 9 . 4 ) ' ) xmax 
call cftext (lxt, lyt, ntex$) 

c 

return 

end 


subroutine mapping (jmax,c) 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
C c 

c Compute the mapping from the grid index I to the grid c 

c locations X. c 

c c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

common / bound / nreg, xb(0:5), rib(0:5), db(0:5), kb(0:5), 

$ rir(0:5), irm(0:5), ril(0:5), ilp(0:5), 

$ xr(0:5), xrm(0:5), xl(0:5), xlp(0:5) 

common / mesh / x(999), dx(999), ddx(999), 

$ imin, imax, itot, xmin, xmax, dxmin, dxmax 


c 

c 

c 


real*8 x(999), dx(999), ddx(999) 
real*8 c(6) 

do 120 n * 1, nreg 

call quintic (n,jmax,c) 

if (kb(n-l) .eq. 1) then 
il ~ ril(n) 
else 

il = ilp(n) 
endif 


if (kb(n) .eq. 1) then 
i2 = rir(n) 
else 

i2 = irm(n) 
endif 
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100 

110 


do 110 i = il, i2 
x(i) = 0 

do 100 j = 2, jmax 
jb = jmax + 2 - j 
x (i) - (x (i) + c ( jb) ) * 

continue 

x ( i ) =x(i) +c(l) 
continue 


120 


continue 


return 

end 


subroutine quintic (n,jmax,c) 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c c 
c Compute the 5th degree polynomial through the end points with c 
c the specified first differences having zero second c 


c 

c 


differences . 


c 

c 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 


common 

/ 

bound 

common 

/ 

mesh 

real*8 

X 

(999) , 

real*8 

a i 

(6,7), 

a (1, 1 ) 

= 

1. 

a (2, 1) 

= 

1. 

a (3, 1) 


1. 

a (4 , 1) 

= 

1. 

do 100 

j 

= 2, j 

a (1, 

j) 

= a (1 


nreg, xb(0:5), rib(0:5), db(0:5), kb(0:5), 
rir (0:5), irm(0:5), ril(0:5), ilp(0:5), 
xr(0:5), xrm(0:5), xl(0:5), xlp(0:5) 
x (999) , dx ( 9 9 9 ) , ddx(999), 

imin, imax, itot, xmin, xmax, dxmin, dxmax 


dx (999) 

s (6, 7) , 


ddx (999) 
c (6) 


100 


a (2, j) 
a (3, j) 
a ( 4 , j) 
continue 


a (2, j-1) 
a (3, j-1) 
a (4, j-1) 


ril(n) 
ilp (n) 
rir (n) 
irm (n) 


200 


do 200 j = 1, jmax 
a (5, j) = ( j - 1) 
a (6, j) = ( j - 1) 
continue 


(j - 2) 
(j - 2) 


ril(n) 
rir (n) 


k k 
k k 


( j 
( j 


- 3) 

- 3) 


a (1,7) 
a (2,1) 
a (3,1) 
a (4,1) 
a (5,1) 
a (6, 7) 


xl (n) 
xlp (n) 
xr (n) 
xrm (n) 

0. 

0. 
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c 

c 


300 

310 

320 

c 

c 

c 


400 

410 


Gaussian elimination Algorithm 6.1 of Burden and Faires 

do 320 jl = 1, jmax - 1 

do 310 j2 = jl + 1, jmax 

s(j2,jl) = a ( j2, jl) / a ( j 1, jl ) 
do 300 j3 = 1, jmax + 1 

a ( j2, j3) = a(j2,j3) - s(j2,jl) * a(jl,j3) 
continue 
continue 
continue 

Backward substitution 

c ( jmax) — a ( jmax, jmax+1) / a (jmax, jmax) 
do 410 jl = 1, jmax - 1 
jb = jmax - jl 
c{ jb) = a ( jb, jmax+1) 
do 400 j2 = jb + 1, jmax 

c ( jb) - c { jb) - a(jb,j2) * c(j2) 
continue 

c( jb) * c( jb) / a( jb, jb) 
continue 

return 

end 


subroutine writegrid (gridfile$) 
c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c ’ c 

c Write the grid to a file in NAMELIST format c 

c c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

common / mesh / x(999), dx(999), ddx(999), 

$ imin, imax, itot, xmin, xmax, dxmin, dxmax 

c 

character ident $*80 
character gridfile$*40 
real*8 x(999), dx(999), ddx{999) 
c 

open (3, file = gridfile$) 

c 

print ' (/ ) ' 

print *, ' You may put an identifying text string on the file.' 

print *, ' Enter the string ' 

read f (a80)', ident$ 

c 

print '(/,a80)', ident$ 

print ' (a) 9 , r GRID = ' 

print r (5 (fll.5, al) ) f , (x(i), i = imin, imax) 

c 

write (3, '(a80) f ) ident$ 
write (3, Ma)M ' GRID - ' 

write (3, ' (5 (f 11 . 5, al) ) ' ) (x(i), i - imin, imax) 

c 

close (3) 

c 

return 

end 
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subroutine outparam 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c c 

c Write the parameters which define the grid to a file. c 

c c 

c The grid is defined by three real numbers. These are c 

c the coordinate, spacing, and grid index value at c 

c each of the grid boundaries. c 

c c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c 


$ 

$ 

$ 


common / bound 


common / mesh 


/ nreg, xb(0:5), rib(0:5), db(0:5), kb(0:5), 
rir (0:5), irm(0:5), ril(0:5), ilp(0:5), 
xr (0 : 5) , xrm(0:5), xl(0:5), xlp(0:5) 

/ x ( 999) , dx ( 9 9 9 ) , ddx(999), 

imin, imax, itot, xmin, xmax, dxmin, dxmax 


real*8 x(999), dx(999), ddx(999) 
open (2, file = 'param.tem') 

write (2, ' (//, a, il, a, $ ) ' ) ' CAPGRID for ', nreg, ' Regions' 

write (2, ' ( 5x, a, i3 ) ' ) ' Itotal = ', itot 

write (2, '(/,3a)') ' Xbound' , ' dXbound' , ' Ibound 

do 100 n = 0, nreg 

write (2, ' (f 11 . 5, f 11 . 5, f 10 . 1) ' ) xb(n), db(n), rib(n) 

100 continue 


close (2) 

return 

end 


subroutine iniparam 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccc 

c c 

c Define default values of the parameters c 

c c 

ccccccccccccccccccccccccccccccccccccccccccccccccccc 


common / bound / nreg, xb(0:5), rib(0:5), db(0:5), kb(0:5), 
$ rir (0:5) , irm(0:5), ril(0:5), i Ip ( 0 : 5 ) , 

$ xr (0 : 5 ) , xrm(0:5), xl(0:5), xlp(0:5) 


if (nreg 

.eq. 1) then 

xb (0 ) 

= 0.00001 

xb (1) 

= 8. 

db (0 ) 

= 0.05 

db (1) 

= 1.00 

rib ( 0 ) 

= 1 . 

rib (1) 

= 20. * 
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elseif (nreg .eq. 2) then 
xb ( 0 ) = 0 . 

xb (1 ) = 1. 

xb (2 ) = 3. 

db (0 ) = 0.18 

db (1) = 0.04 

db (2) = 0.40 

rib (0) = 1.0 

rib(l) = 10.5 
rib (2) = 21.0 
elseif (nreg .eq. 3) then 


xb (0) 

= 

-3. 

xb (1 ) 

- 

0 . 

xb (2 ) 

= 

1 . 

xb (3) 

= 

4. 

db (0) 

= 

0.40 

db (1) 

= 

0.05 

db (2 ) 


0.05 

db (3 ) 

= 

0.40 

rib ( 0 ) 

= 

1.0 

rib (1) 

= 

13.5 

rib (2) 

= 

28.0 

rib (3) 


41 . 0 


elseif (nreg .eq. 4) then 


xb (0 ) 

= 

0.0 

xb (1) 

= 

0.1 

xb (2) 

= 

0.5 

xb (3) 

= 

1.0 

xb (4) 

= 

2.0 

db (0) 

= 

0.04 

db (1) 

= 

0.04 

db (2 ) 

= 

0.03 

db (3) 

= 

0.03 

db (4) 


0.20 

rib (0) 

= 

1.0 

rib (1) 

= 

3.5 

rib (2) 

= 

12 . 5 

rib (3) 

= 

24.5 

rib (4) 


34.0 


elseif (nreg .eq. 5) then 


xb (0) 

= 

-3.0 

xb (1) 


0.0 

xb (2 ) 


1.0 

xb (3) 


1.5 

xb (4 ) 

= 

2.0 

xb (5) 

= 

5.0 

db (0) 


0.50 

db (1) 

= 

0.05 

db (2) 

= 

0.05 

db (3) 

= 

0.05 

db (4) 


0.05 

db (5) 

= 

0.50 

rib ( 0 ) 

= 

1.0 

rib (1) 

= 

11.5 

rib (2) 


24 . 0 

rib (3) 


29.5 

rib (4) 

= 

37.0 

rib (5) 

= 

47.0 


endif 
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c 


return 

end 


subroutine de scrip 
c 

cccccccccccccccccccccccccccccccccccccccecccccc 

c c 

c Print brief instructions on screen c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccc 

c 

character zzz$*8 
c 

print ' {/)' 

print CAP— TSD Grid Design' 

print ' 

print This program will aid in generating a CAP-TSD grid.' 

print The space is divided into separate regions for which' 

print you specify the region boundaries in X units and' 

print give the indices I for these boundaries. These' 

print ' indices may be either integer or half integer' 
print according to whether the boundaries fall on or between' 

print the grid lines, respectively. The program is set up' 

print for a maximum of five regions.' 

print ' 

print You specify the grid spacing at the boundary of each' 

print region. A separate fifth-degree polynomial is used' 

print with zero second derivative at each boundary. Since' 

print the regions are independent you may do a complicated' 

print *, grid one region at a time as long as you use the same' 

print * r boundary values at the common interfaces. ' 

print ' 

print *, The grid and its first and second differences are' 

print listed and plotted. You may write the grid values to' 

print *, ' a file in a NAMELIST format.' 
print *, ' ' 

print The default values of the parameters shown may be used' 
print by responding to the prompt with <Return>. The' 

print * f r program saves your values so that you need not' 
print re-enter parameters that do not change as you improve' 

print your grid. The latest parameters are also saved on' 

print file "param. tern" . ' 

print (/, a, $) ' , ' Enter <Return> to continue ' 
read ' (a8) ' , zzz$ 

return 

end 



OF POOR QLV.UTY 


27 








> rn O Z C 73 Z 

C Z O O z rn 

H O *3 0 

fo m 

(D 00 


